rm(list=ls())
library(tidyverse)
library(modelsummary)
library(jtools)
library(fixest)

library(here)
here::here()  # points to project root


`%nin%` <- Negate(`%in%`)

source("load_data.R")

appropr$pres_party <- as.factor(appropr$pres_party)
appropr$house_party <- as.factor(appropr$house_party)
appropr$sen_party <- as.factor(appropr$sen_party)


panel_grid <- appropr %>%
  mutate(observed = 1) %>%
  complete(comp_acc_id, fiscal_year, fill = list(observed = 0))
panel_grid <- panel_grid %>%
  group_by(comp_acc_id) %>%
  mutate(obs_count = sum(observed)) %>%
  ungroup() %>%
  arrange(obs_count, comp_acc_id) %>%
  mutate(comp_acc_id = factor(comp_acc_id, levels = unique(comp_acc_id)))

panel_structure_plot <- ggplot(panel_grid, aes(x = fiscal_year, y = factor(comp_acc_id))) +
  geom_tile(aes(fill = factor(observed)), color = "white") +
  scale_fill_manual(values = c("0" = "gray90", "1" = "#1B6C92"), labels = c("Not in Appropriations", "In Appropriations")) +
  labs(x = "Fiscal Year", y = "Unit", fill = "Status",
       title = "Panel Structure: science and research accounts in yearly appropriations") +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 0), legend.position = "top")


ggsave("FigS1.pdf", panel_structure_plot, width = 8, height = 4.5)

panel_grid %>% filter(!is.na(department)) %>% dplyr::select(department, program, obs_count, strict) %>% unique() %>% arrange(desc(department), desc(obs_count)) %>% write_csv("TableS12.csv")

appropr %>% group_by(fiscal_year, house_party, sen_party, pres_party) %>% summarise(number_of_accounts = n_distinct(comp_acc_id)) -> yearly_accounts
library(MASS)
mod_num_accounts <- glm.nb(number_of_accounts ~ house_party+ sen_party+ pres_party, yearly_accounts)
summary(mod_num_accounts)
modelsummary(mod_num_accounts, fmt = 2,
             coef_map = c("(Intercept)" = "(Intercept)" ,
                          "house_partyR" = "Rep House",
                          "sen_partyR" = "Rep Senate",
                          "pres_partyR" = "Rep President"), 
             vcov = "HC3",
                          out = "TableS1.docx")



yearly_total_approp <- appropr %>% group_by(fiscal_year) %>% summarise(tot_sci_fund = sum(pl))
yearly_total_approp

toplines <- appropr %>% dplyr::select(fiscal_year, topline, pres_party, house_party, sen_party) %>% unique()

toplines$topline<-toplines$topline*1000

yearly_total_approp <- yearly_total_approp %>% left_join(toplines)

yearly_total_approp <- yearly_total_approp %>%  mutate(share_sci_fund = tot_sci_fund/topline)

sci_share_plot <- ggplot(yearly_total_approp, aes(x=fiscal_year, y=share_sci_fund)) + 
  geom_point(size =2.5) + geom_line(size=2) + theme_minimal(base_size = 15) + 
  geom_col(
           aes(x = fiscal_year, y = .21, fill = pres_party),
           width = 1, alpha = 0.4, show.legend = FALSE) +
  geom_tile(aes(x = fiscal_year, y = .0001, fill = house_party),
            height = 0.13 * diff(range(yearly_total_approp$share_sci_fund)),
            width = 1, show.legend = FALSE) +
  geom_tile(aes(x = fiscal_year, y = .21, fill = sen_party),
            height = 0.13 * diff(range(yearly_total_approp$share_sci_fund)),
            width = 1, show.legend = FALSE) + scale_fill_manual(values = c("#1B6C92","#9B3D26")) +
  # Labels for regions
  annotate("text", x = mean(1980:2020), y = .1,
           label = "Presidency", vjust = 0.5, size = 5, color = "white") +
  annotate("text", x = mean(1980:2020), y = .21,
           label = "Senate", vjust = 0.5, size = 5, color = "white") +
  annotate("text", x = mean(1980:2020), y = 0,
           label = "House", vjust = 0.5, size = 5, color = "white") +
  
  ylab("Science and Research Appropriations\nas Share of Total Budget") + xlab("Fiscal Year")

ggsave("Fig2.pdf", sci_share_plot, width = 8, height = 5)



 ##### Main Appropriation Models

appropr <- appropr %>% mutate(unempl_rate = unempl_rate*100, gdp_per_change = gdp_per_change*100,
                              unempl_rate_lag =  unempl_rate_lag*100, gdp_per_change_lag = gdp_per_change_lag*100)

m_pres_full <- feols(pb/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = appropr)
m_pres_full_controls <- feols(pb/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)


m_hc_full <- feols(hc/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = appropr)
m_hc_full_controls <- feols(hc/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)

m_sc_full <- feols(sc/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = appropr)
m_sc_full_controls <- feols(sc/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)

m_pl <- feols(pl/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = appropr)
m_pl_controls <- feols(pl/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)

m_pl_controls <- feols(pl/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)

vif_out <- car::vif(lm(pl/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr, data = appropr))

write_csv(as_tibble(vif_out, rownames = "Variable"), "TableS12.csv")


##### In Constant dollars
appropr <- appropr %>% mutate(adj_ratio = pb_constant/pb)
appropr <- appropr %>% mutate(topline_constant = topline*adj_ratio)

m_pres_full_constant <- feols(pb_constant/1000000  ~  house_party + sen_party + pres_party + topline_constant | funder, data = appropr)
m_pres_full_controls_constant <- feols(pb_constant/1000000  ~  house_party + sen_party + pres_party + topline_constant  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)


m_hc_full_constant <- feols(hc_constant/1000000  ~  house_party + sen_party + pres_party + topline_constant | funder, data = appropr)
m_hc_full_controls_constant <- feols(hc_constant/1000000  ~  house_party + sen_party + pres_party  + topline_constant + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)

m_sc_full_constant <- feols(sc_constant/1000000  ~  house_party + sen_party + pres_party + topline_constant | funder, data = appropr)
m_sc_full_controls_constant <- feols(sc_constant/1000000  ~  house_party + sen_party + pres_party  + topline_constant + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)

m_pl_constant <- feols(pl_constant/1000000  ~  house_party + sen_party + pres_party  + topline_constant| funder, data = appropr)
m_pl_controls_constant <- feols(pl_constant/1000000  ~  house_party + sen_party + pres_party  + topline_constant + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)

modelsummary(list("Enacted" = m_pl_constant,
                  "Enacted" = m_pl_controls_constant,
                  "House Mark" = m_hc_full_constant,
                  "House Mark" = m_hc_full_controls_constant,
                  "Senate  Mark" = m_sc_full_constant,
                  "Senate  Mark" = m_sc_full_controls_constant,
                  "Pres Request" = m_pres_full_constant,
                  "Pres Request" = m_pres_full_controls_constant), stars = T, 
             coef_map = c("house_partyR" = "Rep House",
                          "sen_partyR" = "Rep Senate",
                          "pres_partyR" = "Rep President",
                          "topline_constant" = "Budget", 
                          "pres_election_year" = "Pres Election Year",
                          "cpi" = "CPI",
                          "unempl_rate" = "Unemployment",
                          "gdp_per_change" = "Δ GDP/pop ",
                          "real_deficit" = "Real Deficit", 
                          "fullyear_cr" = "Full Year CR"), fmt = 2, out = "TableS4.docx")



### Unconditioned on other institutions -- just average differences using OLS

m_pb_pres <- feols(pb/1000000 ~ pres_party  | funder, data = appropr)
m_pb_pres_nofe <- feols(pb/1000000 ~ pres_party ,cluster = ~funder, data = appropr)
m_pb_sen <- feols(pb/1000000 ~ sen_party  | funder, data = appropr)
m_pb_sen_nofe <- feols(pb/1000000 ~ sen_party , cluster = ~funder, data = appropr)
m_pb_house <- feols(pb/1000000 ~ house_party  | funder, data = appropr)
m_pb_house_nofe <- feols(pb/1000000 ~ house_party , cluster = ~funder, data = appropr)

mp_pb <- modelplot(list(m_pb_pres, m_pb_pres_nofe, m_pb_sen, m_pb_sen_nofe, m_pb_house, m_pb_house_nofe), stars = T, 
                 coef_map = c("house_partyR" = "Rep House",
                              "sen_partyR" = "Rep Senate",
                              "pres_partyR" = "Rep President"
                 )) + coord_flip() +scale_color_brewer(palette = "Dark2") + 
  xlab("Estimated Marginal Effect\n(Millions of Dollars)") +
  ggtitle("DV: Presidential Request") + aes(shape = model) +
  scale_shape_manual(values = c(16, 17, 16, 17, 16, 17)) + geom_vline(xintercept = 0, linetype = 3) + theme(legend.position = "none")


m_hc_pres <- feols(hc/1000000 ~ pres_party  | funder, data = appropr)
m_hc_pres_nofe <- feols(hc/1000000 ~ pres_party ,cluster = ~funder, data = appropr)
m_hc_sen <- feols(hc/1000000 ~ sen_party  | funder, data = appropr)
m_hc_sen_nofe <- feols(hc/1000000 ~ sen_party , cluster = ~funder,data = appropr)
m_hc_house <- feols(hc/1000000 ~ house_party  | funder, data = appropr)
m_hc_house_nofe <- feols(hc/1000000 ~ house_party , cluster = ~funder, data = appropr)

mp_hc <- modelplot(list(m_hc_pres, m_hc_pres_nofe, m_hc_sen, m_hc_sen_nofe, m_hc_house, m_hc_house_nofe), stars = T, 
                   coef_map = c("house_partyR" = "Rep House",
                                "sen_partyR" = "Rep Senate",
                                "pres_partyR" = "Rep President"
                   )) + coord_flip() +scale_color_brewer(palette = "Dark2") + 
  xlab("Estimated Marginal Effect\n(Millions of Dollars)") +
  ggtitle("DV: House Mark") + aes(shape = model) +
  scale_shape_manual(values = c(16, 17, 16, 17, 16, 17)) + geom_vline(xintercept = 0, linetype = 3) + theme(legend.position = "none")


m_sc_pres <- feols(sc/1000000 ~ pres_party  | funder, data = appropr)
m_sc_pres_nofe <- feols(sc/1000000 ~ pres_party , cluster = ~funder, data = appropr)
m_sc_sen <- feols(sc/1000000 ~ sen_party  | funder, data = appropr)
m_sc_sen_nofe <- feols(sc/1000000 ~ sen_party , cluster = ~funder, data = appropr)
m_sc_house <- feols(sc/1000000 ~ house_party  | funder, data = appropr)
m_sc_house_nofe <- feols(sc/1000000 ~ house_party , cluster = ~funder, data = appropr)

mp_sc <- modelplot(list(m_sc_pres, m_sc_pres_nofe, m_sc_sen, m_sc_sen_nofe, m_sc_house, m_sc_house_nofe), stars = T, 
                   coef_map = c("house_partyR" = "Rep House",
                                "sen_partyR" = "Rep Senate",
                                "pres_partyR" = "Rep President"
                   )) + coord_flip() +scale_color_brewer(palette = "Dark2") + 
  xlab("Estimated Marginal Effect\n(Millions of Dollars)") +
  ggtitle("DV: Senate Mark") + aes(shape = model) +
  scale_shape_manual(values = c(16, 17, 16, 17, 16, 17)) + geom_vline(xintercept = 0, linetype = 3) + theme(legend.position = "none")


m_pl_pres <- feols(pl/1000000 ~ pres_party  | funder, data = appropr)
m_pl_pres_nofe <- feols(pl/1000000 ~ pres_party , cluster = ~funder, data = appropr)
m_pl_sen <- feols(pl/1000000 ~ sen_party  | funder, data = appropr)
m_pl_sen_nofe <- feols(pl/1000000 ~ sen_party , cluster = ~funder, data = appropr)
m_pl_house <- feols(pl/1000000 ~ house_party  | funder, data = appropr)
m_pl_house_nofe <- feols(pl/1000000 ~ house_party , cluster = ~funder, data = appropr)

mp_pl <- modelplot(list(m_pl_pres, m_pl_pres_nofe, m_pl_sen, m_pl_sen_nofe, m_pl_house, m_pl_house_nofe), stars = T, 
                   coef_map = c("house_partyR" = "Rep House",
                                "sen_partyR" = "Rep Senate",
                                "pres_partyR" = "Rep President"
                   )) + coord_flip() +scale_color_brewer(palette = "Dark2") + 
  xlab("Estimated Marginal Effect\n(Millions of Dollars)") +
  ggtitle("DV: Enacted in Public Law") + aes(shape = model) +
  scale_shape_manual(values = c(16, 17, 16, 17, 16, 17)) + geom_vline(xintercept = 0, linetype = 3) + theme(legend.position = "none")

library(cowplot)
raw_assoc <- plot_grid(mp_pl, mp_pb, mp_hc, mp_sc, labels = "auto")
ggsave("FigS6_nobudgettopline.pdf", width = 7, height = 5)



### Unconditioned on other institutions -- just average differences using OLS -- this time conditioning on the topline

m_pb_pres <- feols(pb/1000000 ~ pres_party + topline | funder, data = appropr)
m_pb_pres_nofe <- feols(pb/1000000 ~ pres_party + topline,cluster = ~funder, data = appropr)
m_pb_sen <- feols(pb/1000000 ~ sen_party + topline | funder, data = appropr)
m_pb_sen_nofe <- feols(pb/1000000 ~ sen_party + topline, cluster = ~funder, data = appropr)
m_pb_house <- feols(pb/1000000 ~ house_party + topline | funder, data = appropr)
m_pb_house_nofe <- feols(pb/1000000 ~ house_party + topline, cluster = ~funder, data = appropr)

mp_pb <- modelplot(list(m_pb_pres, m_pb_pres_nofe, m_pb_sen, m_pb_sen_nofe, m_pb_house, m_pb_house_nofe), stars = T, 
                   coef_map = c("house_partyR" = "Rep House",
                                "sen_partyR" = "Rep Senate",
                                "pres_partyR" = "Rep President"
                   )) + coord_flip() +scale_color_brewer(palette = "Dark2") + 
  xlab("Estimated Marginal Effect\n(Millions of Dollars)") +
  ggtitle("DV: Presidential Request") + aes(shape = model) +
  scale_shape_manual(values = c(16, 17, 16, 17, 16, 17)) + geom_vline(xintercept = 0, linetype = 3) + theme(legend.position = "none")

m_hc_pres <- feols(hc/1000000 ~ pres_party  + topline| funder, data = appropr)
m_hc_pres_nofe <- feols(hc/1000000 ~ pres_party + topline,cluster = ~funder, data = appropr)
m_hc_sen <- feols(hc/1000000 ~ sen_party  + topline| funder, data = appropr)
m_hc_sen_nofe <- feols(hc/1000000 ~ sen_party + topline, cluster = ~funder,data = appropr)
m_hc_house <- feols(hc/1000000 ~ house_party  + topline| funder, data = appropr)
m_hc_house_nofe <- feols(hc/1000000 ~ house_party+ topline , cluster = ~funder, data = appropr)

mp_hc <- modelplot(list(m_hc_pres, m_hc_pres_nofe, m_hc_sen, m_hc_sen_nofe, m_hc_house, m_hc_house_nofe), stars = T, 
                   coef_map = c("house_partyR" = "Rep House",
                                "sen_partyR" = "Rep Senate",
                                "pres_partyR" = "Rep President"
                   )) + coord_flip() +scale_color_brewer(palette = "Dark2") + 
  xlab("Estimated Marginal Effect\n(Millions of Dollars)") +
  ggtitle("DV: House Mark") + aes(shape = model) +
  scale_shape_manual(values = c(16, 17, 16, 17, 16, 17)) + geom_vline(xintercept = 0, linetype = 3) + theme(legend.position = "none")


m_sc_pres <- feols(sc/1000000 ~ pres_party  + topline| funder, data = appropr)
m_sc_pres_nofe <- feols(sc/1000000 ~ pres_party + topline, cluster = ~funder, data = appropr)
m_sc_sen <- feols(sc/1000000 ~ sen_party  + topline| funder, data = appropr)
m_sc_sen_nofe <- feols(sc/1000000 ~ sen_party + topline, cluster = ~funder, data = appropr)
m_sc_house <- feols(sc/1000000 ~ house_party  + topline| funder, data = appropr)
m_sc_house_nofe <- feols(sc/1000000 ~ house_party + topline, cluster = ~funder, data = appropr)

mp_sc <- modelplot(list(m_sc_pres, m_sc_pres_nofe, m_sc_sen, m_sc_sen_nofe, m_sc_house, m_sc_house_nofe), stars = T, 
                   coef_map = c("house_partyR" = "Rep House",
                                "sen_partyR" = "Rep Senate",
                                "pres_partyR" = "Rep President"
                   )) + coord_flip() +scale_color_brewer(palette = "Dark2") + 
  xlab("Estimated Marginal Effect\n(Millions of Dollars)") +
  ggtitle("DV: Senate Mark") + aes(shape = model) +
  scale_shape_manual(values = c(16, 17, 16, 17, 16, 17)) + geom_vline(xintercept = 0, linetype = 3) + theme(legend.position = "none")


m_pl_pres <- feols(pl/1000000 ~ pres_party + topline | funder, data = appropr)
m_pl_pres_nofe <- feols(pl/1000000 ~ pres_party + topline, cluster = ~funder, data = appropr)
m_pl_sen <- feols(pl/1000000 ~ sen_party  + topline| funder, data = appropr)
m_pl_sen_nofe <- feols(pl/1000000 ~ sen_party + topline, cluster = ~funder, data = appropr)
m_pl_house <- feols(pl/1000000 ~ house_party  + topline| funder, data = appropr)
m_pl_house_nofe <- feols(pl/1000000 ~ house_party + topline, cluster = ~funder, data = appropr)

mp_pl <- modelplot(list(m_pl_pres, m_pl_pres_nofe, m_pl_sen, m_pl_sen_nofe, m_pl_house, m_pl_house_nofe), stars = T, 
                   coef_map = c("house_partyR" = "Rep House",
                                "sen_partyR" = "Rep Senate",
                                "pres_partyR" = "Rep President"
                   )) + coord_flip() +scale_color_brewer(palette = "Dark2") + 
  xlab("Estimated Marginal Effect\n(Millions of Dollars)") +
  ggtitle("DV: Enacted in Public Law") + aes(shape = model) +
  scale_shape_manual(values = c(16, 17, 16, 17, 16, 17)) + geom_vline(xintercept = 0, linetype = 3) + theme(legend.position = "none")


raw_assoc <- plot_grid(mp_pl, mp_pb, mp_hc, mp_sc, labels = "auto")
ggsave("FigS6.pdf", width = 7, height = 5)


#### As share of the topline
appropr <- appropr %>% mutate(pl_share = pl/topline, pb_share = pb/topline, hc_share = hc/topline, sc_share = sc/topline)


m_pres_full_share <- feols(pb_share  ~  house_party + sen_party + pres_party  | funder, data = appropr)
m_pres_full_controls_share <- feols(pb_share  ~  house_party + sen_party + pres_party   + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)


m_hc_full_share <- feols(hc_share  ~  house_party + sen_party + pres_party  | funder, data = appropr)
m_hc_full_controls_share <- feols(hc_share  ~  house_party + sen_party + pres_party   + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)

m_sc_full_share <- feols(sc_share  ~  house_party + sen_party + pres_party  | funder, data = appropr)
m_sc_full_controls_share <- feols(sc_share  ~  house_party + sen_party + pres_party   + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)

m_pl_share <- feols(pl_share  ~  house_party + sen_party + pres_party  | funder, data = appropr)
m_pl_controls_share <- feols(pl_share  ~  house_party + sen_party + pres_party   + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)

modelsummary(list("Enacted" = m_pl_share,
                  "Enacted" = m_pl_controls_share,
                  "House Mark" = m_hc_full_share,
                  "House Mark" = m_hc_full_controls_share,
                  "Senate  Mark" = m_sc_full_share,
                  "Senate  Mark" = m_sc_full_controls_share,
                  "Pres Request" = m_pres_full_share,
                  "Pres Request" = m_pres_full_controls_share), stars = T, 
             coef_map = c("house_partyR" = "Rep House",
                          "sen_partyR" = "Rep Senate",
                          "pres_partyR" = "Rep President",
                          "topline" = "Budget", 
                          "pres_election_year" = "Pres Election Year",
                          "cpi" = "CPI",
                          "unempl_rate" = "Unemployment",
                          "gdp_per_change" = "Δ GDP/pop ",
                          "real_deficit" = "Real Deficit", 
                          "fullyear_cr" = "Full Year CR"), fmt = 2, out = "TableS5.docx")




#### Aggregated to the agency level
appropr %>% group_by(agency_id, fiscal_year) %>% summarise(pb_agency = sum(pb),
                                                           hc_agency = sum(hc),
                                                           sc_agency = sum(sc),
                                                           pl_agency = sum(pl)) -> agency_totals 

appropr %>% dplyr::select(agency_id, fiscal_year, house_party ,sen_party , pres_party , topline  , pres_election_year , cpi , unempl_rate , gdp_per_change , real_deficit , fullyear_cr) %>% unique() -> covars 
agency_totals <- agency_totals %>% left_join(covars)
agency_totals

m_pres_full_agency <- feols(pb_agency/1000000 ~  house_party + sen_party + pres_party  + topline | agency_id, data = agency_totals)
m_pres_full_controls_agency <- feols(pb_agency/1000000  ~  house_party + sen_party + pres_party  + topline + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| agency_id, data = agency_totals)


m_hc_full_agency <- feols(hc_agency/1000000  ~  house_party + sen_party + pres_party + topline | agency_id, data = agency_totals)
m_hc_full_controls_agency <- feols(hc_agency/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| agency_id, data = agency_totals)

m_sc_full_agency <- feols(sc_agency/1000000  ~  house_party + sen_party + pres_party  + topline| agency_id, data = agency_totals)
m_sc_full_controls_agency <- feols(sc_agency/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| agency_id, data = agency_totals)

m_pl_agency <- feols(pl_agency/1000000  ~  house_party + sen_party + pres_party + topline | agency_id, data = agency_totals)
m_pl_controls_agency <- feols(pl_agency/1000000  ~  house_party + sen_party + pres_party  + topline + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| agency_id, data = agency_totals)


modelsummary(list("Enacted" = m_pl_agency,
                  "Enacted" = m_pl_controls_agency,
                  "House Mark" = m_hc_full_agency,
                  "House Mark" = m_hc_full_controls_agency,
                  "Senate  Mark" = m_sc_full_agency,
                  "Senate  Mark" = m_sc_full_controls_agency,
                  "Pres Request" = m_pres_full_agency,
                  "Pres Request" = m_pres_full_controls_agency), stars = T, 
             coef_map = c("house_partyR" = "Rep House",
                          "sen_partyR" = "Rep Senate",
                          "pres_partyR" = "Rep President",
                          "topline" = "Budget", 
                          "pres_election_year" = "Pres Election Year",
                          "cpi" = "CPI",
                          "unempl_rate" = "Unemployment",
                          "gdp_per_change" = "Δ GDP/pop ",
                          "real_deficit" = "Real Deficit", 
                          "fullyear_cr" = "Full Year CR"), fmt = 2, out = "TableS6.docx")



#### With size of the majority
appropr$house_rep_margin <- appropr$house_rep - appropr$house_dem
appropr$sen_rep_margin <- appropr$sen_rep - appropr$sen_dem

m_pres_margin <- feols(pb/1000000  ~  house_rep_margin + sen_rep_margin + pres_party + topline | funder, data = appropr)
m_pres_margin_controls <- feols(pb/1000000  ~  house_rep_margin + sen_rep_margin + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)


m_hc_margin <- feols(hc/1000000  ~  house_rep_margin + sen_rep_margin + pres_party + topline | funder, data = appropr)
m_hc_margin_controls <- feols(hc/1000000  ~  house_rep_margin + sen_rep_margin + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)

m_sc_margin <- feols(sc/1000000  ~  house_rep_margin + sen_rep_margin + pres_party + topline | funder, data = appropr)
m_sc_margin_controls <- feols(sc/1000000  ~  house_rep_margin + sen_rep_margin + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)

m_pl_margin <- feols(pl/1000000  ~  house_rep_margin + sen_rep_margin + pres_party + topline | funder, data = appropr)
m_pl_margin_controls <- feols(pl/1000000  ~  house_rep_margin + sen_rep_margin + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = appropr)

modelsummary(list("Enacted" = m_pl_margin,
                  "Enacted" = m_pl_margin_controls,
                  "House Mark" = m_hc_margin,
                  "House Mark" = m_hc_margin_controls,
                  "Senate  Mark" = m_sc_margin,
                  "Senate  Mark" = m_sc_margin_controls,
                  "Pres Request" = m_pres_margin,
                  "Pres Request" = m_pres_margin_controls), stars = T, 
             coef_map = c("house_rep_margin" = "Rep House Margin",
                          "sen_rep_margin" = "Rep Senate Margin",
                          "pres_partyR" = "Rep President",
                          "topline" = "Budget", 
                          "pres_election_year" = "Pres Election Year",
                          "cpi" = "CPI",
                          "unempl_rate" = "Unemployment",
                          "gdp_per_change" = "Δ GDP/pop ",
                          "real_deficit" = "Real Deficit", 
                          "fullyear_cr" = "Full Year CR"), fmt = 2, output = "TableS8.docx")

######### Only on "strict" definition of science

m_pres_full_strict <- feols(pb/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = filter(appropr, strict == 1))
m_pres_full_strict_controls <- feols(pb/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = filter(appropr, strict == 1))


m_hc_full_strict <- feols(hc/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = filter(appropr, strict == 1))
m_hc_full_strict_controls <- feols(hc/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = filter(appropr, strict == 1))

m_sc_full_strict <- feols(sc/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = filter(appropr, strict == 1))
m_sc_full_strict_controls <- feols(sc/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = filter(appropr, strict == 1))

m_pl_strict <- feols(pl/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = filter(appropr, strict == 1))
m_pl_strict_controls <- feols(pl/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = filter(appropr, strict == 1))

modelsummary(list("Enacted" = m_pl_strict,
                  "Enacted" = m_pl_strict_controls,
                  "House Mark" = m_hc_full_strict,
                  "House Mark" = m_hc_full_strict_controls,
                  "Senate  Mark" = m_sc_full_strict,
                  "Senate  Mark" = m_sc_full_strict_controls,
                  "Pres Request" = m_pres_full_strict,
                  "Pres Request" = m_pres_full_strict_controls), stars = T, 
             coef_map = c("house_partyR" = "Rep House",
                          "sen_partyR" = "Rep Senate",
                          "pres_partyR" = "Rep President",
                          "topline" = "Budget", 
                          "pres_election_year" = "Pres Election Year",
                          "cpi" = "CPI",
                          "unempl_rate" = "Unemployment",
                          "gdp_per_change" = "Δ GDP/pop ",
                          "real_deficit" = "Real Deficit", 
                          "fullyear_cr" = "Full Year CR"), fmt = 2, out = "TableS7.docx")




#### Funding main model

m_funding <- feols(funding_usd/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = filter(funding, start_year > 1979))

#Effect sizes aggregated 
coef(m_funding)["house_partyR"]*m_funding$fixef_sizes*1000000

coef(m_funding)["sen_partyR"]*m_funding$fixef_sizes*1000000

coef(m_funding)["pres_partyR"]*m_funding$fixef_sizes*1000000

library(equatiomatic)
summary(m_funding)
#### main models table
modelsummary(list("Enacted" = m_pl,
                  "Enacted" = m_pl_controls,
                  "House Committee Mark" = m_hc_full,
                  "House Committee Mark" = m_hc_full_controls,
                  "Senate Committee Mark" = m_sc_full,
                  "Senate Committee Mark" = m_sc_full_controls,
                  "Pres Request" = m_pres_full,
                  "Pres Request" = m_pres_full_controls,
                  "Grantmaking" = m_funding), stars = T, 
             coef_map = c("house_partyR" = "Rep House",
                          "sen_partyR" = "Rep Senate",
                          "pres_partyR" = "Rep President",
                          "topline" = "Budget", 
                          "pres_election_year" = "Pres Election Year",
                          "cpi" = "CPI",
                          "unempl_rate" = "Unemployment",
                          "gdp_per_change" = "Δ GDP/pop ",
                          "real_deficit" = "Real Deficit", 
                          "fullyear_cr" = "Full Year CR"), fmt = 2, out = "TableS2.docx")


mp1 <- modelplot(list("Enacted" = m_pl,
                      "Enacted Appropriations with Controls" = m_pl_controls),
             coef_map = c("house_partyR" = "Rep House",
                          "sen_partyR" = "Rep Senate",
                          "pres_partyR" = "Rep President"
                          )) + coord_flip() +
  scale_color_brewer("", palette = "Dark2") + 
  xlab("Estimated Marginal Effect\n(Millions of Dollars)") + geom_vline(xintercept = 0, linetype = 3) +theme(legend.position = "bottom")


mp2 <- modelplot(list(
               "Full Institutional Model" = m_hc_full,
               "Full Institutional Model with Controls" = m_hc_full_controls),
          coef_map = c("house_partyR" = "Rep House",
                       "sen_partyR" = "Rep Senate",
                       "pres_partyR" = "Rep President"
          )) + coord_flip() +scale_color_brewer(palette = "Dark2") + xlab("Estimated Marginal Effect\n(Millions of Dollars)") +
  ggtitle("DV: House Committee Mark") + geom_vline(xintercept = 0, linetype = 3) +theme(legend.position = "bottom")


mp3 <- modelplot(list( 
                "Full Institutional Model" = m_sc_full,
                "Full Institutional Model with Controls" = m_sc_full_controls), stars = T, 
          coef_map = c("house_partyR" = "Rep House",
                       "sen_partyR" = "Rep Senate",
                       "pres_partyR" = "Rep President"
          )) + coord_flip() +scale_color_brewer(palette = "Dark2") + 
  xlab("Estimated Marginal Effect\n(Millions of Dollars)") +
  ggtitle("DV: Senate Committee Mark")+ geom_vline(xintercept = 0, linetype = 3) +theme(legend.position = "bottom")


mp4 <- modelplot(list(
               "Full Institutional Model" = m_pres_full, 
               "Full Institutional Model with Controls" = m_pres_full_controls), stars = T, 
          coef_map = c("house_partyR" = "Rep House",
                       "sen_partyR" = "Rep Senate",
                       "pres_partyR" = "Rep President"
          )) + coord_flip() +scale_color_brewer(palette = "Dark2") + 
  xlab("Estimated Marginal Effect\n(Millions of Dollars)") +
  ggtitle("DV: Presidential Request") + geom_vline(xintercept = 0, linetype = 3) +theme(legend.position = "bottom")

mp_funding <- modelplot( m_funding, 
  coef_map = c("house_partyR" = "Rep House",
               "sen_partyR" = "Rep Senate",
               "pres_partyR" = "Rep President"
  )) + coord_flip() +scale_color_brewer(palette = "Dark2") + 
  xlab("Estimated Marginal Effect\n(Millions of Dollars)") +
  ggtitle("Outgoing Grants") + geom_vline(xintercept = 0, linetype = 3) +theme(legend.position = "bottom")


ggsave("FigS5A.pdf", mp_funding, device = "pdf", width = 4, height = 4, units = "in")

budgetstage_plot <- plot_grid(mp4, mp2, mp3, labels = "AUTO", nrow = 1)
ggsave("FigS4.pdf", budgetstage_plot, device = "pdf", width = 11, height = 4, units = "in")

ggsave("Fig3.pdf", mp1, device = "pdf", width = 11, height = 4, units = "in")

####### Descriptive plots


funding_sum <- funding %>% 
  group_by(fiscal_year) %>%
  summarise(num_funders = length(unique(funder)), 
            total_funding = sum(funding_usd_adj))

appropr_sum <- appropr %>% 
  group_by(fiscal_year) %>% 
  summarise(num_accounts = length(unique(funder)), 
            total_approps = sum(pl_constant))


appropr_avg <- appropr %>% group_by(program, department) %>% summarise(total_pl = sum(pl_constant), num_years = length(unique(fiscal_year)), avg_appropr = total_pl/num_years)

dept_appropr <- appropr %>% group_by(department, fiscal_year) %>% summarise(total_pl = sum(pl_constant)) %>% ungroup() %>% group_by(fiscal_year) %>% mutate(fund_rank = dense_rank(desc(total_pl)))
dept_appropr$department_label <- dept_appropr$department

dept_appropr$department_label[dept_appropr$fund_rank > 7] <- "Other"

dept_appropr <- dept_appropr %>% group_by(department_label, fiscal_year) %>% summarise(total_pl = sum(total_pl))


library(directlabels)
p_line <- ggplot(dept_appropr, aes(x=fiscal_year, y=total_pl, group=department_label, color = department_label, label = department_label)) + geom_line() + theme_minimal() + scale_color_brewer(palette = "Dark2") + xlim(1980,2030) + scale_y_continuous("Yearly Appropriations", labels = scales::dollar)

p_line <-  direct.label(p_line, method = list("last.points"))


library(treemapify)

appropr_topdept <- appropr %>% group_by(department, fiscal_year) %>% summarise(total_dept_pl = sum(pl_constant))
appropr_topdept <- appropr_topdept %>% ungroup() %>% group_by(fiscal_year) %>% mutate(fund_rank = dense_rank(desc(total_dept_pl)))

appropr<- appropr %>% left_join(appropr_topdept)
appropr$department_label <- appropr$department

appropr$department_label[appropr$fund_rank > 7] <- "Other"


tm20 <- ggplot(filter(appropr, fiscal_year == 2020), aes(area = pl_constant, fill = department_label, subgroup = department_label, label = program)) +
  geom_treemap() + theme(legend.position = "none") + geom_treemap_subgroup_border() +
  geom_treemap_subgroup_text(place = "centre", grow = T, alpha = 0.5, colour =
                               "black", fontface = "italic", min.size = 0) + scale_fill_brewer(palette = "Dark2")




library(cowplot)


FigS2BC <- plot_grid(p_line, tm20, labels = c("B", "C"), nrow = 1)
ggsave("FigS2BC.pdf", FigS2BC, device = "pdf", width = 10, height = 7)







##### Run permutation test script now

source("permutations.R")

bar_pvals <- readRDS("adj_pvals.rds")
bar_pvals


party_group_funding <- appropr %>% group_by(house_party, major_group) %>% summarise(num_years = length(unique(fiscal_year)),
                                                                                    total_funding = sum(pl_constant, na.rm = T),
                                                                                    yearly_funding = total_funding/num_years)  



party_group_funding$major_group[is.na(party_group_funding$major_group)] <- "Other"

party_group_funding$major_group <- fct_reorder(party_group_funding$major_group, party_group_funding$yearly_funding) %>% fct_rev()

party_group_funding$adj_pvals <- NA
party_group_funding$adj_pvals[1:10] <- bar_pvals



party_group_funding <- party_group_funding %>% mutate(sig_symbol = case_when(adj_pvals >.10 ~ "ns",
                                                                             adj_pvals <=.001 ~ "***",
                                                                             adj_pvals <=.01 ~ "**",
                                                                             adj_pvals <=.05 ~ "*",
                                                                             adj_pvals <=.10 ~ "."))


party_funding_bar <- party_group_funding  %>% ggplot(aes(y=major_group, x=yearly_funding/1000000, fill = house_party)) + 
  geom_bar(stat="identity", position = "dodge") + geom_text(aes(label =sig_symbol, x = -3000)) +
  scale_fill_manual("Party Control of the House", values = c("#156B90", "#9A3E25")) +
  theme(legend.position = "bottom") + xlab("Average Appropriations per Fiscal Year (1980-2020)\n(in Millions of 2021 Constant Dollars)") + ylab("") +
  scale_x_continuous(labels=scales::dollar) + theme_minimal() + theme(legend.position = "top", axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=1)) +
  coord_flip()

ggsave("Fig4.pdf", party_funding_bar, device = "pdf", width = 6, height = 6)




#### Estimated separately by major agency groups
appropr$major_group[is.na(appropr$major_group)] <- "Other"


approp_pl_subsets <- appropr %>% split(.$major_group) %>%
  map(~feols(pl/1000000 ~  house_party + sen_party + pres_party + topline | funder, data = .)) 

approp_pl_subsets_dat <- approp_pl_subsets %>% 
  map(~tidy(.x)) %>%
  bind_rows(.id = "FundingGroup") %>% 
  filter(term!="topline") %>%
  mutate(
         term = case_when(term == "house_partyR" ~ "Republican Control of House",
                          term == "pres_partyR" ~ "Republican President",
                          term == "sen_partyR" ~ "Republican Control of Senate"), 
         lower.outer = estimate - 1.96*std.error,
         upper.outer = estimate + 1.96*std.error,
         lower.inner = estimate - 1.44*std.error,
         upper.inner = estimate + 1.44*std.error)

approp_pl_subsets_dat$FundingGroup <- fct_reorder(approp_pl_subsets_dat$FundingGroup, approp_pl_subsets_dat$estimate, "max") %>% fct_rev()


modelsummary(approp_pl_subsets,
             stars = T, 
             coef_map = c("house_partyR" = "Rep House",
                          "sen_partyR" = "Rep Senate",
                          "pres_partyR" = "Rep President",
                          "topline" = "Budget"),
             output = "TableS3.docx")


Departmental_Approps_Plot <- ggplot(approp_pl_subsets_dat, aes(y=estimate,
                              x=FundingGroup)) + 
  geom_point() + 
  geom_errorbar(aes(ymin = lower.inner, ymax = upper.inner), width = 0, position = position_dodge(width=.5), size = 2, alpha=.3) +
  geom_errorbar(aes( ymin = lower.outer, ymax = upper.outer), width = 0, position = position_dodge(width=.5), size = 1, alpha=.3) + 
  theme_minimal()+ scale_y_continuous(labels = scales::dollar) + 
  scale_color_brewer(type = "qual", palette = "Set2") + 
  geom_hline(yintercept = 0, linetype = 3) +
  facet_wrap(~term, ncol = 1) + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=1), legend.position = "none") +
  xlab("")  +ylab("Estimated Marginal Effect")
  
ggsave("FigS3.pdf", Departmental_Approps_Plot, device = "pdf", width = 6, height = 6, units = "in")




m_pres_full_pre2000 <- feols(pb/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = filter(appropr, fiscal_year <= 2000))
m_pres_full_controls_pre2000  <- feols(pb/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = filter(appropr, fiscal_year<=2000))


m_hc_full_pre2000  <- feols(hc/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = filter(appropr, fiscal_year<=2000))
m_hc_full_controls_pre2000  <- feols(hc/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = filter(appropr, fiscal_year<=2000))

m_sc_full_pre2000  <- feols(sc/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = filter(appropr, fiscal_year<=2000))
m_sc_full_controls_pre2000  <- feols(sc/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = filter(appropr, fiscal_year<=2000))

m_pl_pre2000  <- feols(pl/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = filter(appropr, fiscal_year<=2000))
m_pl_controls_pre2000  <- feols(pl/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = filter(appropr, fiscal_year<=2000))


modelsummary(list("Enacted" = m_pl_pre2000 ,
                  "Enacted" = m_pl_controls_pre2000 ,
                  "House Committee Mark" = m_hc_full_pre2000 ,
                  "House Committee Mark" = m_hc_full_controls_pre2000 ,
                  "Senate Committee Mark" = m_sc_full_pre2000 ,
                  "Senate Committee Mark" = m_sc_full_controls_pre2000 ,
                  "Pres Request" = m_pres_full_pre2000 ,
                  "Pres Request" = m_pres_full_controls_pre2000 ),
             coef_map = c("house_partyR" = "Rep House",
                          "sen_partyR" = "Rep Senate",
                          "pres_partyR" = "Rep President",
                          "topline" = "Budget", 
                          "pres_election_year" = "Pres Election Year",
                          "cpi" = "CPI",
                          "unempl_rate" = "Unemployment",
                          "gdp_per_change" = "Δ GDP/pop ",
                          "real_deficit" = "Real Deficit", 
                          "fullyear_cr" = "Full Year CR"), fmt = 2, stars = T, output = "TableS9.docx")


m_pres_full_post2000 <- feols(pb/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = filter(appropr, fiscal_year > 2000))
m_pres_full_controls_post2000  <- feols(pb/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = filter(appropr, fiscal_year > 2000))


m_hc_full_post2000  <- feols(hc/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = filter(appropr, fiscal_year >2000))
m_hc_full_controls_post2000  <- feols(hc/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = filter(appropr, fiscal_year > 2000))

m_sc_full_post2000  <- feols(sc/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = filter(appropr, fiscal_year > 2000))
m_sc_full_controls_post2000  <- feols(sc/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = filter(appropr, fiscal_year > 2000))

m_pl_post2000  <- feols(pl/1000000  ~  house_party + sen_party + pres_party + topline | funder, data = filter(appropr, fiscal_year>2000))
m_pl_controls_post2000  <- feols(pl/1000000  ~  house_party + sen_party + pres_party + topline  + pres_election_year + cpi + unempl_rate + gdp_per_change + real_deficit + fullyear_cr| funder, data = filter(appropr, fiscal_year > 2000))


modelsummary(list("Enacted" = m_pl_post2000 ,
                  "Enacted" = m_pl_controls_post2000 ,
                  "House Committee Mark" = m_hc_full_post2000 ,
                  "House Committee Mark" = m_hc_full_controls_post2000 ,
                  "Senate Committee Mark" = m_sc_full_post2000 ,
                  "Senate Committee Mark" = m_sc_full_controls_post2000 ,
                  "Pres Request" = m_pres_full_post2000 ,
                  "Pres Request" = m_pres_full_controls_post2000 ),
             coef_map = c("house_partyR" = "Rep House",
                          "sen_partyR" = "Rep Senate",
                          "pres_partyR" = "Rep President",
                          "topline" = "Budget", 
                          "pres_election_year" = "Pres Election Year",
                          "cpi" = "CPI",
                          "unempl_rate" = "Unemployment",
                          "gdp_per_change" = "Δ GDP/pop ",
                          "real_deficit" = "Real Deficit", 
                          "fullyear_cr" = "Full Year CR"), fmt = 2, stars = T, output = "TableS10.docx")

